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LARGE  EDDY  SIMULATIONS  OF  COMPLEX  TURBULENT 
FLOWS  USING  IMMERSED  BOUNDARY  METHOD 


MAYANK  TYAGI,  SUMANTA  ACHARYA 
Mechanical  Engineering  Department 
Louisiana  State  University,  Baton  Rouge,  LA 

Abstract 

Large  eddy  simulations  (LES)  are  performed  for  two  representative 
complex  geometry  problems  of  particular  interest  to  turbomachinery  flows.  First 
problem  studied  is  the  flow  field  inside  a trapped-vortex  combustor.  Second 
problem  is  the  unsteady  interaction  of  rotor  blade  with  stator  blade  wake  field. 
The  representation  of  complex  geometry  is  done  using  immersed  boundary 
method  (IBM).  Two  different  implementations  are  presented  for  the  body  force 
terms. 

1.  Introduction 

Simulation  of  turbulent  flows  in  complex  geometries  is  a daunting  task. 
LES  can  formally  alleviate  the  issue  of  ever-increasing  resolution  demand  for 
high  Reynolds  number  flow.  However,  complex  geometries  pose  the  problem  of 
commutation  errors  on  curvilinear  grids.  Moreover,  the  representation  of  moving 
geometries  using  either  sliding  meshes  or  regenerating  the  mesh  becomes 
overwhelmingly  complicated  in  complex  situations.  IBM  relies  upon  the  body 
force  terms  added  in  the  momentum  equations  to  represent  the  geometry  on  a 
fixed  Cartesian  mesh  (Peskin,  1977,  Mohd.-Yusof,  1996,  Glowinski  et  ah,  1994, 
Fadlun  et  al.  2000,  Kellog,  2000).  This  formulation  is  simple  and  ideally  suited 
for  the  moving  geometries  involving  no-slip  walls  with  prescribed  trajectories 
and  locations. 


2.  Governing  Equations 

In  the  immersed  boundary  method,  the  complex  geometrical  features  are 
incorporated  by  adding  a forcing  function  in  the  governing  equations.  The 
forcing  function  is  zero  everywhere  except  at  the  surface  where  the  influence  of 
the  solid  boundaries  is  assigned  (Subscript  T). 
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where  the  convective  terms  are  represented  by  C and  the  diffusion  terms  are 
represented  by  D.  In  case  of  highly  refined  meshes,  it  may  be  necessary  to  treat 
some  directions  implicitly  for  diffusion  terms  (generally  using  Crank-Nicholson 
scheme). 


C=-(m- V)m,Z)=— V2m 
Re 

To  obtain  the  pressure  Poisson  equation,  take  the  divergence  of  the  second  step 
and  enforce  the  continuity  condition  for  the  velocity  field  at  the  next  time  step 
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Therefore,  the  Poisson  equation  for  pressure  can  be  solved  prior  to  the  second 
step  in  the  time-split  scheme.  The  spatial  discretization  is  performed  using  fourth 
order  central  difference  schemes.  For  the  two  problems  studied,  two  different 
implementations  for  the  body  force  terms  are  presented.  The  trapped-vortex 
combustor  is  a complicated  but  stationary  geometry  and  hence,  interpolation  on 
only  one-side  of  the  geometry  is  adequate.  However,  for  moving  rotor  blade  the 
interpolation  is  performed  on  both  sides  of  the  curved  surface. 


CASE  A:  Forcing  at  only  one  side  of  the  immersed  boundary  (inside  the  virtual 
solid) 
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Figure  1 Identification  of  the  circular  boundary  on  uniform  2-D  cartesian  mesh 
and  evaluation  of  the  nearest  exterior  point  corresponding  to  each  identified 
interior  point 

Let  A be  the  mesh  spacing  and  S be  the  distance  of  the  forcing  point  from  the 
immersed  surface.  Therefore,  we  apply  the  linear  interpolation/extrapolation 
among  the  forced  point,  point  on  the  immersed  surface  and  the  point  just  outside 
the  virtual  solid.  Let  Vd  be  the  desired  velocity  at  the  point  on  the  immersed 


LES  OF  COMPLEX  TURBULENT  FLOWS  USING  IBM  41 9 

surface  and  Vc  be  the  computed  velocity  in  the  region  of  interest.  Therefore,  the 
velocity  at  the  forcing  point  Vim  is  given  by 

[Vc  - Vd]/[  A - 8]  = [Vd  - Vim]/[5] 

Vim  = Vd  [A]/[  A - 8]  - Vc  [8]/[  A - 8] 

Clearly,  In  the  limit  8 going  to  zero,  i.e.  the  forcing  point  approaching  the  point 
on  the  immersed  surface,  we  retrieve  the  limit  Vim  approaching  Vd.  However,  In 
the  limit  8 approaching  mesh  spacing  A,  we  have  Vc  approaching  Vd.  Vjm  is  ill- 
defined  because  it  is  the  difference  between  Vd  and  Vc  with  very  large 
coefficients.  For  such  coefficients,  Vc  is  set  equal  to  Vd. 

CASE  B:  Forcing  at  both  sides  of  the  immersed  boundary  (inside  the  virtual 
solid  and  at  the  very  first  point  outside  the  virtual  solid) 

Let  A be  the  mesh  spacing  and  8 be  the  distance  of  the  forcing  point  from 
the  immersed  surface.  Therefore,  we  apply  the  linear  interpolation/extrapolation 
among  the  forced  points,  point  on  the  immersed  surface  and  the  computed  point 
just  outside  the  virtual  solid  (Figure  2).  Let  Vd  be  the  desired  velocity  at  the  point 
on  the  immersed  surface  and  Vc  be  the  computed  velocity  in  the  region  of 
interest.  Therefore,  the  velocity  at  the  internal  forcing  point  VlM  is  given  by 

[Vc-Vd]/[2A-S]  = [Vd-Vint]/[8] 

Vim  = Vd  [2A]/[  2A  - 8]  - Vc  [8]/[  2A  - 8] 

Similarily,  the  velocity  at  the  external  forcing  point  Vext  is  given  by 

[Vc  - Vd]/[  2A  - 8]  = [Vex,  - Vd]/[A  - 8] 

Vext  = Vd  [A]/[  2A  - 8]  + Vc  [A  - 8]]/[  2A  - 8] 

Clearly,  in  the  limit  8 going  to  zero,  i.e.  the  internal  forcing  point 
approaching  the  point  on  the  immersed  surface,  we  retrieve  the  limit  Vim 
approaching  Vd  and  Vexi  approaching  [Vd  + Vc]/2.  Moreover,  In  the  limit  8 
approaching  mesh  spacing  A,  we  have  Vext  approaching  Vd . V;nt  is  defined  as 
[2Vd  - Vc]  as  it  should  be  by  reflection  condition.  Thus,  the  forcing  remains 
physical  for  all  positions  of  the  immersed  surface  between  the  grid  interfaces. 
For  moving  boundary  implementation,  it  will  be  important  for  another  reason.  It 
avoids  reflected  velocities  to  show  up  in  the  region  of  interest. 
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Figure  2 Identification  of  the  interpolation  stencils  for  a moving  rotor  blade. 
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3.  Trapped  Vortex  Combustor 

3.1  Problem  Description 

A uniform  cartesian  grid  of  92x57x117  points  is  used  for  a domain 
containing  the  upper  half  of  the  TV  combustor.  All  the  dimensions  are  selected  to 
approximate  the  experimental  setup  of  Mancilla  (2001)  within  Cartesian 
geometry.  Ratio  of  air  injection  velocity  to  the  mainflow  velocity  is  2.2. 
Reynolds  number  based  on  the  annular  mainflow  velocity  and  air  hole  dimension 
(D)  is  3400  for  these  simulations.  The  radii  of  the  forebody,  the  connecting  tube, 
the  afterbody  and  the  outer  shell  are  24.5D,  3.7D,  23D  and  27.5D  respectively. 
The  lengths  of  the  forebody,  the  connecting  tube  and  the  afterbody  are  1 2D,  30D 
and  12D  respectively.  The  periodicity  of  the  geometry  is  exploited  by  putting 
half  of  the  jet  injections  around  the  bottom  of  the  computational  domain  with  the 
boundary  conditions  obtained  by  the  rotational  symmetry  about  the  axial 
direction.  At  the  inflow,  fully  developed  laminar  profile  along  with  fluctuations 
is  prescribed.  The  fluctuations  are  assumed  to  be  Gaussian  and  are  calculated 
using  Box-Muller  algorithm.  At  the  walls,  no  slip  boundary  conditions  are 
imposed  using  immersed  boundary  method.  Uniform  injections  of  air  and  fuel 
are  applied  at  the  respective  holes  on  the  afterbody.  Periodic  boundary  conditions 
are  applied  in  the  spanwise  (z)  direction.  At  the  outflow,  a non-reflective 
convective  scheme  is  applied  to  convect  away  the  flow  structures  out  of  the 
computational  domain  without  any  spurious  reflections.  The  wave  speed  is 
calculated  to  maintain  the  mass  flux  balance  in  the  whole  domain. 

3.2  Results 

The  three-dimensional  simulations  showed  that  the  vorticity  magnitude 
of  the  trapped  vortex  in  the  cavity  changes  due  to  vortex  stretching  mechanism 
that  is  absent  in  the  two-dimensional  simulations  done  by  other  researchers 
(Katta  and  Roquomore,  1996  and  Stone  and  Menon,  2000).  The  motion  of  the 
trapped  vortex  is  unsteady  inside  the  cavity  as  observed  experimentally  as  well  as 
numerically  (Mancilla,  2001).  Unsteady  dynamics  of  the  coherent  structures 
inside  the  cavity  is  a slower  process  as  compared  to  the  separation  region  over 
the  afterbody  and  the  mixing  layer  region  behind  the  forebody  lip.  The  ingestion 
of  annular  mainflow  in  front  of  the  afterbody  separation  region  is  the  main 
mechanism  of  the  flow  entrainment  inside  the  cavity.  The  mixture  in  the  cavity  is 
ejected  radially  outwards  due  to  the  pressure  gradients  there.  The  instantaneous 
axial  and  radial  locations  of  the  TV  are  different  at  different  meridional  planes. 
However,  time-averaged  streamtrace  patterns  seem  to  converge  towards  a TV 
with  simpler  geometrical  features. 

TIME  AVERAGED  FLOW  FIELD 

The  streamtraces  are  presented  at  the  meridional  planes  0 = 90°,  0°  and 
180°  (figure  3).  At  0 = 90°,  the  large  recirculation  region  is  formed  between  the 
annular  mainflow  and  the  fuel  injections.  The  unsteadiness  of  the  mixing  layer 
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and  the  motion  of  TV  inside  the  cavity  has  been  averaged  out.  At  0 = 0°  and 
180°,  the  streamtraces  are  asymmetric  implying  that  a longer  averaging  period  is 
needed  to  completely  eliminate  the  slow  variations.  Moreover,  the  axisymmetry 
is  also  absent  at  these  meridional  planes.  The  TV  is  a doughnut  shaped  structure 
inside  the  cavity.  The  pitch  of  the  spiral  formed  by  streamtraces  close  to  the 
meridional  plane  0 = 90°  is  largest.  The  azimuthal  motion  of  the  streamtraces 
along  the  surface  of  the  TV  is  absent  in  all  the  previous  2-D  studies.  However, 
the  core  of  the  TV  is  mostly  irrotational. 


Time  averaged  alreamllnes 


(b) 

Figure  3 Streamtraces  at  (a)  0 = 90°  (b)  0 = 0°  and  180°  (tracers  are 
released  along  radial  positions  at  several  axial  locations). 

TURBULENT  STRESSES  (No  Figures) 

At  0 = 90°,  the  turbulent  shear  stress  u’v’  inside  the  cavity  is  mostly 
generated  by  the  production  terms  ( u’u’.dV/dx  and  v V.dU/dy)  around  the 
mixing  layer  region.  The  high  axial  momentum  fluctuations  of  the  fluid  parcels 
are  correlated  with  radially  inward  fluctuations  around  the  mixing  layer  region. 
This  implies  the  turbulent  mixing  and  diffusion  is  enhanced  between  the  cavity 
and  annular  mainflow  by  this  stress  component.  The  shear  stress  is  also  large 
along  the  jet  injections.  However,  the  levels  of  turbulent  stresses  indicate  that 
flow  is  mostly  laminar  close  to  the  forebody  and  connecting  tube  junction.  At  0 = 
0°  and  180°,  the  turbulent  shear  stress  u’w’  inside  the  cavity  is  mostly  generated 
by  the  production  terms  (u’u’.dW/dx  and  w’w’’.3U/3z)  around  the  jets  in  the 
cavity  and  the  mixing  layer  region.  This  distribution  indicates  enhanced  mixing 
of  fluid  parcels  inside  the  core  of  TV.  The  distribution  of  u’w’  near  the  jet 
injections  is  dictates  the  radial  spread  of  the  jets  along  the  axial  direction  in  the 
cavity.  At  X/D  = 40.8,  the  radial  fluctuations  are  high  at  different  azimuthal  or 
meridional  planes.  The  normal  stress  u’u’  is  mostly  negligible  due  to  suppression 
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of  fluctuations  normal  to  the  face  of  the  afterbody.  The  only  significant  levels  are 
above  the  afterbody  in  the  annular  region.  The  distribution  of  the  shear  stress 
v’w’  along  the  jets  is  primarily  due  to  turbulent  production  (v’v’.3W/3y  and 
w'w’  .dV/dz).  The  distribution  in  the  annular  region  seems  to  attain  the 
periodicity  between  the  fuel  injections.  It  can  be  expected  that  the  modeling  of 
the  shear  stress  tensor  components  using  an  eddy  viscosity  approximation  may 
work  here  since  the  shear  stress  components  depend  on  the  corresponding  mean 
strain  rate  tensor  components  only. 

4.  Unsteady  Stator-Rotor  Interactions 

4.1  Problem  Description 

The  inherent  unsteadiness  of  a turbomachinery  flow  field  created  by 
relative  motion  between  stationary  blades  (stator)  and  the  rotating  blades  (rotor), 
requires  the  designer  to  account  for  three-dimensional  as  well  as  unsteady 
effects.  The  unsteadiness  is  caused  by  (a)  the  interaction  of  the  rotor  airfoils  with 
the  wakes  and  passage  vortices  generated  by  upstream  airfoils,  (b)  the  relative 
motion  of  the  rotors  with  respect  to  the  stators  (potential  effect),  and  (c)  the 
shedding  of  vortices  by  the  airfoils  because  of  the  blunt  trailing  edges  (Rai  and 
Madavan,  1990,  Saxer  and  Giles,  1994).  Computation  of  such  flows  is 
complicated  by  relative  motion  between  rotor  and  stator  airfoils  and  the  periodic 
transition  of  the  flow  from  laminar  to  turbulent.  Unsteady  simulations  have  been 
performed  using  multitudes  of  approximations  such  as  unsteady  RANS, 
“average-passage”  approach  and  “mixing-plane”  approach.  These  calculations 
have  been  performed  invariably  using  “sliding  mesh”  techniques  requiring 
further  modeling  of  “apparent  stresses”  and  “phase-lagged”  interface  conditions 
(Sharma  et  al.,  1994).  In  the  present  study,  we  utilize  LES  with  moving  IBM  to 
simulate  unsteady  stator-rotor  interactions.  Though,  the  calculation  is  performed 
for  incompressible  fluid  at  a low  Reynolds  number,  it  demonstrates  the  strength 
of  the  method  by  avoiding  all  ad-hoc  assumptions  pertaining  to  RANS  modeling 
and  sliding  meshes.  A uniform  Cartesian  grid  of  302x202x1 1 points  is  used  for  a 
domain  of  the  size  3DxlDx0.1D,  where  D is  the  chord  length  of  the  rotor  airfoil. 
Choice  of  a small  spanwise  dimension  may  not  allow  larger  physical  scales  and 
hence  may  not  be  desirable.  The  geometry  of  the  airfoils  is  taken  from  the 
numerical  study  of  Kelecy  et  al  (1995).  The  airfoil  profile  is  approximated  by  the 
cubic  spline  surfaces.  The  airfoil  is  divided  into  leading  edge,  trailing  edge, 
pressure  surface  and  suction  surface  to  ensure  that  immersed  boundary 
conditions  are  enforced  on  enough  grid  points  to  realize  the  geometry.  Uniform 
flow  field  is  specified  at  the  inflow.  Periodic  boundary  conditions  are  applied  in 
the  direction  of  rotor  motion  (y)  and  the  spanwise  (z)  direction.  Reynolds  number 
base  on  the  inflow  velocity  and  rotor  chord  length  is  5000.  At  the  outflow,  a non- 
reflective  convective  scheme  is  applied  to  convect  away  the  flow  structures  out 
of  the  computational  domain  without  any  spurious  reflections.  The  wave  speed  is 
calculated  to  maintain  the  mass  flux  balance  in  the  whole  domain. 
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4.2  Results 

Four  snapshots  of  the  instantaneous  vorticity  field  are  shown  in  figure  4 
at  time  instants  T apart  (=0.1 2x,  where  x is  the  non-dimensional  time  scale). 
There  is  a separation  region  on  the  suction  side  of  the  stator.  The  trailing  edge 
vortices  of  the  stator  blade  impact  on  the  suction  side  of  the  rotor  blade  near  its 
leading  edge.  The  trailing  edge  vortices  of  the  rotor  and  the  vortices  formed  due 
to  the  interaction  of  stator  wake  and  suction-side  boundary  layer  are  shed  into  the 
passage  flow  and  convected  out  of  the  domain. 


d)  to  + 3T 


Figure  4 Instantaneous  snapshots  of  vorticity  component  at  time  a)  t = to,  b)  t = 
t0+T,  c)  t = t0+2T  and  d)  t = t0+3T. 
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5.  Concluding  Remarks 

The  merit  of  immersed  boundary  method  (IBM)  for  simulating  moving 
complex  geometries  on  Cartesian  mesh  has  been  demonstrated.  High  order  of 
accuracy  of  discretization  schemes  is  retained  which  is  very  important  for  LES. 
The  problems  studied  are  a)  Trapped  vortex  combustor  (TVC)  and  b)  Stator-rotor 
interactions  in  a transonic  turbine  stage.  In  TVC  case  study,  mixing  inside  an 
annular  cavity  is  analyzed.  In  stator-rotor  case  study,  the  superiority  of  this 
method  is  demonstrated  over  existing  methodologies  such  as  sliding  meshes. 
Moreover,  the  ad-hoc  modeling  for  the  “^parent  stresses”  is  not  needed  in  the 
realms  of  LES.  In  future,  a zonal  refinement  treatment  of  the  immersed 
boundaries  will  be  implemented  to  capture  the  essential  near  wall  physics  to 
render  this  method  with  predictive  capabilities. 
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